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Abstract 

Multistep matrix splitting iterations serve as preconditioning for Krylov subspace meth¬ 
ods for solving singular linear systems. The preconditioner is applied to the generalized min¬ 
imal residual (GMRES) method and the flexible GMRES (FGMRES) method. We present 
theoretical and practical justifications for using this approach. Numerical experiments show 
that the multistep generalized shifted splitting (GSS) and Hermitian and skew-Hermitian 
splitting (HSS) iteration preconditioning are more robust and efficient compared to standard 
preconditioners for some test problems of large sparse singular linear systems. 
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1 Introduction 

Consider solving linear systems 


Ax = 6, (1.1) 

where A € M nxn may be singular and b € M n . For solving large sparse linear systems (1.1), 
iterative methods are preferred to direct methods in terms of efficiency and memory require¬ 
ment. When the problem (1.1) is ill-conditioned, the convergence of iterative methods such 
as Krylov subspace methods tends to deteriorate and may be accelerated by using precondi¬ 
tioning. However, well-established preconditioners using incomplete matrix factorizations [33], 
[4], [9] require additional memory whose amount is typically comparable to that of the given 
problem, and may not work in the singular case. 

Another approach for preconditioning Krylov subspace methods for solving linear systems is 
to use a splitting matrix such as the successive overrelaxation (SOR) method [22], [46]. Matrix 
splitting iterations can serve as preconditioning for Krylov subspace methods. 

In the singular case, some iterative methods and preconditioners may be infeasible, i.e., 
they may break down and/or fail to converge. In this paper, we focus on using GMRES with 
preconditioning since the method is well-established and fairly well understood in the singular 
case [49], [23], [21]. GMRES applied to the linear system (1.1) with initial iterate x$ € gives 
the fcth iterate x such that ||6 — Aa^H = min 3 . Ga . 0 _|_£ fe (A,r 0 ) ||6 —A*||, where || • || is the Euclidean 
norm, = b — Ax o is the initial residual, and JCk(A, r*o) = spanjro, Atq, ..., A k ~ l ro] is the 
Krylov subspace of order k. Hereafter, denote /C& = ICk(A,r q) for simplicity. 
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In the singular case, GMRES may fail to determine a solution of (1.1). GMRES is said 
to break down at some step k if dimA/Cfc < dim /Cfc or dim/Cfc < k [11, p. 38]. Note that 
dimA/Cfc < dim/Cfc < k holds for each k. The dimensions of A/Cfc and /Cfc are related to the 
uniqueness of the iterate £Cfc, whereas dim/Cfc is related to the degeneracy of the Krylov subspace 
method. GMRES determines a solution of Ax = b without breakdown for all b € 7C(A) and 
for all xo € M n if and only if A is a group (GP) matrix A 7(A) fl 7 Z(A) = {0} [11, Theorem 2.6], 
[30, Theorem 2.2], cf. [35, Theorem 4.4.6], where A 7(A) is the null space of A and 71(A) is the 
range space of A. The condition that A is a GP matrix is equivalent to that the largest size of 
the Jordan block of A corresponding to eigenvalue 0 is not larger than one [31, section 3]. 

Other than Krylov subspace methods, much efforts have been made to study matrix split¬ 
ting iterations for solving singular linear systems (1.1) (see [27], [29], [18], [10], [36], [47], [37], 
[15], [41]). Some of modern matrix splitting iterations were shown to be effectively used as pre¬ 
conditioning for Krylov subspace methods, and can be potentially useful as multistep matrix 
splitting iteration preconditioning. For example, see [3], [28], [17], [44], [45] for the Hermitian 
and skew-Hermitian splitting (HSS) iterations, [43] for the triangular and skew symmetric split¬ 
ting (TSS) iterations, [14] for the generalized shift splitting (GSS) iterations, and [50], [48] for 
Uzawa methods for singular saddle point problems. We shed some light on the preconditioning 
aspect of matrix splitting iterations in the singular case. 

Consider applying GMRES to the preconditioned linear system AP~ 1 u = b, x = P" 1 u, 
which is equivalent to Ax = b, where P is nonsingular and a preconditioning matrix given 
by multistep matrix splitting iterations. The right-preconditioned GMRES (RP-GMRES) 
method with initial iterate xq € M n determines the A:th iterate x fc such that ||b — Aa?fc|| = 
min a . ga . 0+ y Cfc (p-i J 4 j p-i ro ) ||6 — Aic||, where uq € M n and r$ = b — AP~ 1 uq = b — Ax q. On the 
other hand, the flexible GMRES (FGMRES) method [34] allows to change the precondition¬ 
ing matrix for each iteration. This means that the number of the multistep matrix splitting 
iterations may vary in GMRES. 

The rest of the paper is organized as follows. In section 2, we give sufficient conditions 
such that GMRES preconditioned by a fixed number of matrix splitting iterations determines 
a solution without breakdown, a spectral analysis of the preconditioned matrix, and a con¬ 
vergence bound of the method, and discuss the computational complexity of the method. In 
section 3, we give sufficient conditions such that FGMRES preconditioned by multistep matrix 
splitting iterations determines a solution without breakdowns in the singular case. In section 
4, we show numerical experiment results on test problems comparing the multistep generalized 
shift-splitting (GSS) and Hermitian and skew-Hermitian (HSS) matrix splitting iteration pre¬ 
conditioners with the GSS and HSS preconditioners, respectively. In section 5, we conclude the 
paper. 

2 GMRES preconditioned by a fixed number of matrix splitting 
iterations 

Consider applying a preconditioner using several steps of matrix splitting iterations to RP- 
GMRES. We give its algorithm as follows (cf. [19], [20]). 

Here, C (£) is the preconditioning matrix given by a fixed number i of matrix splitting 
iterations, e* is the zth column of the identity matrix, and 77 m +i.m = {hi,j} € 

We next give an expression for the preconditioned matrix AC^ for GMRES with i matrix 
splitting iterations. Consider the matrix splitting iterations applied to Az = in line 3. Note 
Ufc € 71(A) if b G 71(A). Let M be a nonsingular matrix such that A = M—N. Denote the itera¬ 
tion matrix by H = M~ 1 N. Assume that the initial iterate is € J\f(H), e.g., z ^ = 0. Then, 
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Algorithm 2.1 GMRES method preconditioned by I matrix splitting iterations. 

1: Let xq £ R” be the initial iterate, tq :=b — Ax o, /3 := ||r*o||, Vi '■= tq/ 

2: for k = 1,2,... until convergence do 

3: Apply I iterations of a matrix splitting to Az = v k to obtain z k = C^v k ', 

4: w := Az k , for i = 1,2, ... ,k do hi )k ; = w := w — hi t kVi end for 

5: if hfc+i k := ||m|| = 0 then set m := k and go to line 7 else v k+ \ := w/h k +\ k', 

6: end for 

7: y m ■■= argmin yeRm ||/3ei - H m+1:m y\\, x m := x Q + [z 1 ,z 2 ,.. .,z m \y m ; 


the lih iterate of the matrix splitting iterations is z^ = Hz ^ ^ + M l v k = E;=o RiM lv k, 
I £ N. Hence, the multistep matrix splitting iteration preconditioning and preconditioned 
matrices are C® = E*L ( 1 L/W -1 and AC® = M~ l E^d H 'i l ~ H ) M = M ~\ l ~ Rt ) M = 
M" 1 Ei=o(I — M~ l A) 1 AI~ X AM , respectively. 

We will give sufficient conditions such that GMRES preconditioned by matrix splitting 
iterations determines a solution of (1.1) without breakdown. First, we prepare the following 

Proposition 2.1 ([24], [31, Theorem 1], [38, Theorem 2]). Let H be a square real matrix. Then, 
H is semiconvergent, i.e., lim*-^ H l exists, if and only if either A = 1 is semisimple, i.e., the 
algebraic and geometric multiplicities corresponding to A = 1 are equal, or |A| < 1 holds for all 
A € cr(H) = {A | Hv = Xv,v € C n \{0}} the spectrum of H. 

Lemma 2.2. If H is semiconvergent, then Ef=o nonsingular for all I € N. 

Proof. Proposition 2.1 shows that there exists a nonsingular matrix S such that J = S~ l HS = 
J © I is the Jordan canonical form (JCF) of H with p(J) < 1 for J € M rxr , where © denotes 
the direct sum and p(A) = max{|A| : A £ cr(A)} is the spectral radius of A. Hence, E;Lo ^ = 
S {[(I — J) -1 (I - J 1 )] © (H)| 5 _1 holds for all i € N. Since 1 - A £ / 0 holds for all A £ a(J) 
and for all I £ N, I — is nonsingular and hence Ei=o is nonsingular for all i £ N. □ □ 

Lemma 2.3. If H is semiconvergent, then I — H { is a GP matrix for all I £ N. 

Proof. If O is the zero matrix, then I — H ( = S'[(I — .T ) © 0]5 _1 . Since I — is nonsingular, 
I — H l is a GP matrix for all l £ N. □ □ 

Now we show that GMRES preconditioned by a fixed number of matrix splitting iterations 
determines a solution of Ax = b. 

Theorem 2.4. Assume that the iteration matrix H is semiconvergent. Then, GMRES pre¬ 
conditioned by multistep matrix splitting iterations defined above determines a solution of 
Ax = b without breakdown for all b £ 'JZ(A), for all xq £ M n , and for all I £ N. 

Proof. Since Ei=o nonsingular for all I € N from Lemma 2.2, C^ = E;=o H l M~ 1 is 

nonsingular for all I £ N. Hence, the preconditioned linear system C^Ax = C^b is equivalent 
to Ax = b. Since C^A = I — H f is a GP matrix for all I £ N from Lemma 2.3, the theorem 
follows from [30, Theorem 2.2], □ □ 

This theorem gives [30, Theorem 4.6] as a corollary if the linear system (1.1) is a symmetric 
and positive semidefinite linear system. 

Theorem 2.4 relies on the property that the preconditioned matrix is GP, which is implied 
by the semiconvergence of the iteration matrix, irrespective of the property of A. Hence, we 
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may extend the class of singular linear systems that GMRES can solve by combining with pre¬ 
conditioners. This means that even though A is not a GP matrix, the multistep matrix splitting 
iteration preconditioned matrix is a GP matrix for H semiconvergent, and GMRES precondi¬ 
tioned by the multistep matrix splitting iterations determines a solution without breakdown 
(see Theorem 5.1 in Appendix). 

Semiconvergence is a simple and convenient property for deciding if a matrix splitting 
method is feasible as multistep matrix splitting iterations for preconditioning GMRES in the 
singular case. Indeed, there are many matrix splitting iterations whose iteration matrix can be 
semiconvergent. They are powerful when used as matrix splitting preconditioners for Krylov 
subspace methods, and potentially useful as multistep matrix splitting iteration preconditioning 
for GMRES such as the Jacobi, Gauss-Seidel, SOR, and symmetric SOR (SSOR) methods [18], 
extrapolated methods [36], two-stage methods [41], the GSS method [14], and the HSS method 
and and its variants [3], [28], [17], [44], [45]. 

Theorem 2.4 applies to some trivial examples. For example, if A = L + D + L J is symmetric 
and positive semidefinite and the SOR splitting matrix is M = + wL), where D is 

diagonal, L is strictly lower triangular, and u: € R, then the SOR iteration matrix H = M _1 N 
is semiconvergent for cj € (0,2) [18, Theorem 13]. On the other hand, if an iteration matrix 
H is semiconvergent, the extrapolated iteration matrix (1 — y)I + jH is also semiconvergent 
for 0 < 7 < 2/(1 + v(H)) [36, Theorem 2.2], We will recall conditions such that the GSS and 
HSS iteration matrices are semiconvergent in sections 4.1 and 4.2, respectively. Hence, these 
multistep matrix splitting iterations can serve as preconditioning for GMRES. 


2.1 Spectral analysis and convergence bound 


Next, consider decomposing GMRES preconditioned by multistep matrix splitting iterations 
into the IZ(AC^) = 1Z(A) and IZ(A) 1 - components to analyze the spectral property of the pre¬ 
conditioned matrix (cf. [23]). Assume that the iteration matrix H is semiconvergent through¬ 
out this subsection and b € 1Z{A). Let r = rankA, Q\ € R nxr SU ch that 1Z(Qi) = 'JZ(A), 
Q 2 € R nx ( n_r ) such that 1Z(Q2) = TZ(A)- 1 , and Q = [Qi,Q 2 ] be orthogonal. Then, GMRES 
applied to AC^u = b can be seen as GMRES applied to (Q T AC^Q)Q T u = Q J b , or 


'QjAC^Qi 

qJacWq 2 ' 

Qju 


[An 

A 12 

u 1 


'QJb 


V 

0 

0 

Qju 


0 

0 

u 2 


Qjb 


0 


Assume that the initial iterate satisfies uq € 7Z(A). Then, the kth iterate of GMRES applied 
to ( 2 . 1 ) is given by 


Q 'Hk 



eQ J u 0 + Q J IC k (AC^,r 0 ) 


Ur 


+ 


An A12 

O O 



which minimizes ||Q T (fe — AC^u k ) |, or uj, e Uq + /C k (An, rj) which minimizes ||rfc|| = H5 1 — 
Anu \||. This means that u), is equal to the kth iterate of GMRES applied to Anu 1 = b 1 for 
all k (cf. [23, section 2.5 ]). 

Now we give the spectrum of the preconditioned matrix AC^ to present a convergence 
bound on GMRES preconditioned by multistep matrix splitting iterations. Let r = rank A. 
The r nonzero eigenvalues of AC^ are the eigenvalues of An, since 


det ( AC w - Al) 


det 


( An — AI r 

U O 


A12 
A I n—r 


(—X) n ~ r det(An - AI r ) 


and An is nonsingular AC^' is a GP matrix [23, Theorem 2.3]. If yi is an eigenvalue of 
H , then AC^ = M^ 1 (I — H l )M has an eigenvalue A = 1 — //. From Proposition 2.1, H has 
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r eigenvalues such that |/i| < 1 and n — r eigenvalues such that y = 1. For \p\ < 1, we obtain 
|A — 11 = \p.\ e < u(H) e < 1, where n(H) = max{|A| : A € cr(iL)\{l}} is the pseudo spectral 
radius of H, i.e., the r eigenvalues of AC ^ are in a disk with center at 1 and radius u(HY < 1. 
For p = 1, we have A = 0, i.e., the remaining n — r eigenvalues are zero. 

Theorem 2.5. Let r k be the kth residual of GMRES preconditioned by t matrix splitting it¬ 
erations and T be the Jordan basis of AC^\ Assume that H is semi-convergent. Then, 
we have ||rfc|| < k(T) El=o’^ Ci) P(-^) ki ~ l \\ r o\\ f or x o € 7 Z(C^A) and for all b € 11(A), 
where k(T) = ||T|| ||T _1 ||, d is the size of the largest Jordan block corresponding to a nonzero 
eigenvalue of C^A, and r(k,d) = min (k,d— 1). 

Proof. Theorem 2.4 ensures that GMRES preconditioned by multistep matrix splitting itera¬ 
tions determines a solution of Ax = b without breakdown for all b € 1Z(A) and for all xq € R n . 
From [1, Theorem 1], we have 

\\ r k\\ = min \\p(AC^) r o\\ < k(T) min max |b(^)llll r o||, 

pePfc,p(o)=i pePfc,p(o)=i i<i<s 


where P*. is the set of all polynomials of degree not exceeding k and J* is a Jordan block of AC 1 ' 1 ' 1 
corresponding to a nonzero eigenvalue, i = 1,2,..., s. The second factor is bounded above by 
min peP fe ,p(o)=i maxi<j< s \\p(Ji)\\ < Eilo’ d) ( k )p( H ) U ~ l [ , Theorems 2, 5]. □ □ 

Note that the residual |||| does not necessarily depend only on the eigenvalues of AC^ 
when k(T) is large (see [39] and references therein). 

2.2 Computational complexity 

Compare GMRES for k£ iterations preconditioned by one step of a matrix splitting method 
with that for k iterations preconditioned by l matrix splitting iterations of the same matrix 
splitting method in terms of the Krylov subspaces for the iterate x k e, x k and computational 
complexity, since they use the same total number of matrix splitting iterations. 

Proposition 2.6. If C ^ and H are as defined above and H is semiconvergent, then we have 

K k (C^A,C^r 0 ) C IC ke (C^A,C^r 0 ). 

Proof. The proof is by induction. Let A = M~ ] A and fo = Consider the case k = 1. 

We have /Ci(C^l4, C^Vo) = spanjC^Vo} and 

C®r 0 = £(I - M-'AfM-'ro = £ £ ( l \ (-i) J f 0 . 

i =0 *=0 j =0 V - V 

Since C^tq € spanjfo, Ar o,..., A e ~ 1 f , o}, we have 

/Ci(C (£ U,C w r 0 ) C /Q(C ( 1 l4,C ( 1 ) r 0 ). 

Next, assume that K, k (C^A,C^ro) C tC k e(C^A, C^tq) = JCke(A,ro) holds. Then, 
JC k+1 (CWA, C®r 0 ) = KL k (C^A, C®r Q ) + span {(C^A) k C^r 0 }, 

%+i )t(C (1) A, cWro) = K m (A, f 0 ) + IC k (A, A kl r 0 ). 
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From C^A = Y^\= o(I — v4)*^4, we have 


(C^A) k C w r 0 


n-i 


k +1 


£(I - Aj 


A k r 0 , 


which belongs to lC(k+i)e(A,ro) = JC^ k+ i\£(C^A,C^ro). Hence, 

K k+1 {C®A, C®r 0 ) C JC ke (CWA, C«r 0 ) = /C (fc+1) ^(i, f 0 ) holds. □ □ 

This proposition shows that GMRES preconditioned by one matrix splitting iteration gives 
an optimal Krylov subspace for the iterate, i.e., any Krylov subspace given by GMRES for k 
iterations preconditioned by t matrix splitting iterations is not larger than the one given by 
GMRES for kl iterations preconditioned by one matrix splitting iteration. However, GMRES 
preconditioned by one matrix splitting iteration is not necessarily more efficient than GMRES 
preconditioned by more than one matrix splitting iteration, as will be seen in section 4. In¬ 
deed, while GMRES for kl iterations preconditioned by one matrix splitting iteration requires 
kt matrix-vector products of A with z k and kl orthogonalizations, GMRES for k iterations 
preconditioned by l matrix splitting iterations requires k matrix-vector products of A with z k 
and k orthogonalizations. Hence, GMRES for kl iterations preconditioned by one matrix split¬ 
ting iteration needs more computations. Therefore, GMRES preconditioned by more than one 
matrix splitting iteration may be more efficient. 

Moreover, Proposition 2.6 gives a lower bound of the number of iterations of GMRES precon¬ 
ditioned by multistep matrix splitting iterations which is required to determine a solution. Let 
s be the smallest integer such that ]C S (C^A, C^ro) < s. Assume that GMRES preconditioned 
by l matrix splitting iterations determines a solution at the kth step, where k is the smallest 
integer such that K. k (C^A,C^ro) = K. s -i(C^A, C^ro). Then, k is larger than (s — 1 )/£. 
Hence, GMRES preconditioned by £ matrix splitting iterations requires more than (s — 1 )/£ 
iterations to determine a solution. 

3 Flexible GMRES preconditioned by multistep matrix split¬ 
ting iterations. 

The preconditioners given in section 2 uses a fixed number of matrix splitting iterations. This 
can be extended to allow a variable number of matrix splitting iterations for each iteration in line 
3, Algorithm 2.1 (flexible GMRES (FGMRES) method [34]). Let be the multistep matrix 
splitting iteration preconditioning matrix for the kill iteration. Then, the FGMRES iterate x\ 
is determined over the space xq + TZ(Z^) = xq + ,..., C^Av\]), which is 

no longer a Krylov subspace. Quantities denoted with superscript F are relevant to FGMRES 
hereafter. Hence, Theorem 2.4 does not apply to FGMRES preconditioned by multistep matrix 
splitting iterations. 

Similarly to the breakdown of GMRES due to the linear dependence of v k+ \ on uj, « 2 , ■ • •, 
v k , FGMRES may break down with /if +1 k = 0 due to the matrix-vector product AC^ k = °> 
i.e., v\ € A f(AC^A), in the singular case. If is nonsingular, then for b e TZ(A), / 0 
-€=> AC^Av ¥ k / 0 is equivalent to that AC^A is a GP matrix, which is given by the iteration 
matrix H semiconvergent. 

Notice that [34, Proposition 2.2] holds irrespective of the nonsingularity of A: if tq ^ 0, 
hj +l i / 0 for * = 1, 2 ,... , k — 1, and = {hf j} € U. kxk is nonsingular, then h\ +lk = 0 
is equivalent to that the FGMRES iterate x\ is uniquely determined and is a solution of 
Ax = b. Here, the nonsingularity of H* is ensured by an additional assumption as follows. 
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Let Qk J Rk+i,k = fc be the QR factorization of H k+l k , where Qk is the product of Givens 
rotations ■ ■ ■ IR such as 1R = R_i © *•] © R-* and Rk+i,k £ R( fc +i)xfc is upper 

triangular. The scalars Ck and Sk are chosen to satisfy c k + s k = 1 and to vanish the ( k + 1, k) 
entry of hR_i • • • fR H k+1 k . It follows from [40, Lemma 4] that if ||uj^ — Az \|| < |c^— 1 1 for c\ 0, 
C 2 0, ..., Ck- 1 0 and r\ 0, then Hk is nonsingular. Thus, we have the following. 

Theorem 3.1. If the iteration matrix H defined above is semiconvergent and the multistep 
matrix splitting iterations attain the residual norm ||r;^ —Az^|| < |ca,| for the kth iteration, then 
FGMRES preconditioned by the multistep matrix splitting iterations determines a solution of 
Ax = b for all b € R(A) and for all Xq € R n . 


4 Numerical experiments 

Numerical experiments on the discretized Stokes problem and artificially generated problems 
show the feasibility of GMRES and FGMRES preconditioned by multistep matrix splitting 
iterations and the effectiveness of the former. These methods were compared with previous 
preconditioners in terms of the central processing unit (CPU) time. For instance for multi- 
step matrix splitting iteration preconditioning, we used the generalized shift-splitting (GSS) 
and Hermitian and skew-Herrnitian splitting (HSS) and their inexact variants. Although no 
condition such that GMRES preconditioned by a fixed number of matrix splitting iterations of 
an inexact splitting determines a solution without breakdown is given, we used the method for 
comparisons. 

The initial iterates for the multistep matrix splitting iterations and GMRES and FGMRES 
iterations were set to zero. No restarts were used for these methods. The matrix splitting 
iterations in FGMRES approximately solved the linear system Az = Vk to the accuracy on 
the residual norm — Az]! +1 || < |q.| to ensure that FGMRES determines a solution without 
breakdown (Theorem 3.1). The stopping criterion used for GMRES and FGMRES iterations 
was in terms of the relative residual norm ||b — Accfc|| < 10 _6 ||ro||. 

The computations were done on a computer with an Intel Xeon CPU E5-2670 2.50GHz, 256 
GB random-access memory (RAM), and Community Enterprise Operating System (CentOS) 
Version 6.8. All programs for the iterative methods were coded and run in Matlab R2014b for 
double precision floating point arithmetic with unit roundoff 2~ 53 ~ 1.1 ■ 10~ 16 . 


4.1 Multistep generalized shifted splitting iteration preconditioning. 

We give numerical experiment results on singular saddle point problems 


Ax = 


C B T 
-B O 


x = b, Be R qxp , C e R pxp positive definite, 


(4.1) 


comparing GMRES preconditioned by I iterations and FGMRES preconditioned by £k iterations 
of the generalized shifted splitting (GSS) 


1 

0:1 + C 

B J 

z (i+ 1) — I 

cd 

-B t 

2 

—B 


2 

B 

fil _ 


W+d, i= 1,2,or 4 (4.2) 


and its inexact variant (IGSS) with GMRES with the standard GSS and IGSS preconditioning 
£ = 1 and a sparse direct solver, where £ is the number of GSS and IGSS iterations. 

Consider test problems of the form (4.1) given by the Stokes problem — /j,Au + Vp = /, 
V • u = 0 in an open domain Q in R 2 with the boundary and normalization conditions u = 0 
on and f n p(x)dx = 0, respectively, where p is the kinematic viscosity constant, A is the 


7 








componentwise Laplace operator, the vector field u denotes the velocity, V and V- denote the 
gradient and divergence operators, respectively, and the scalar function p denotes the pressure. 
The Stokes problem was discretized upwind in square domain 17 = (0,1) x (0,1) on uniform 
grid. Thus, the matrix representation of the Stokes problem is C = (I 9 < 8 > T + T <g> I 9 ) © (I g <S> T + 
T®Iq) € R 2 <? 2 x 2 g 2 ; £T = [BT )6ljb2 ] G R 2q 2 x(g 2 +2)^ = [(I g <g, i?) T , (F ® I ? ) T ] € W?* 2 ? ,T = 
/x/i _ 2 tridiag(—1, 2, —l) + (2/i) _ 1 tridiag(—1,1,0) € M. qxg ,F = /i _ 1 tridiag(—1,1,0) € W xg , where 
<g> denotes the Kronecker product, bj = [e T ,0 T ]R, bj = [0 T ,e T ]R, e = [1,1, ..., 1] T € IF 2 / 2 , 
and h = (q + l)” 1 is the discretization meshsize [8]. The (2,1) and (1,2) blocks were modified 
to be rank-deficient as done in [50, section 5], [14, Example 4.1]. We chose two viscosity values 
p, = 10 -5 and 1 and three kinds of grids 16 x 16, 24 x 24, and 32 x 32. The right-hand side 
vector for ( 1 . 1 ) was set to b = Ae. 

The GSS iteration matrix is semiconvergent for a, f3 > 0 [14, Theorem 3.2] and GMRES 
preconditioned by the multistep GSS iterations determines a solution of Ax = b without break¬ 
down for (4.1), since A is positive definite (Theorem 2.4). On the other hand, the GMRES 
methods preconditioned by IGSS and its multistep version are not guaranteed to determine a 
solution without breakdown. The value of f3 for GSS and IGSS was set to ||R|| 2 /||C'|| [13]. The 
value of a for GSS and IGSS was experimentally determined to have the minimal CPU time. 
The resulting values were a = 10, 13, and 15 for /x = 1 with grids 16, 24, and 32, respectively, 
and a = 30, 37, and 57 for fx = 10 -5 with grids 16, 24, and 32, respectively. 

The linear system (4.2) was solved via [13, Algorithm 2.1] by using the LU factorization for 
GSS and was solved by using GMRES with the stopping criterion 10 _1 in terms of the relative 
residual norm for the inexact multistep GSS (IGSS) iteration preconditioning [7, Section 6 ]. 

Tables 4.1 and 4.2 give the number of iterations and the CPU time in seconds for the Stokes 
problem with /i = 1 and 10~ 5 , respectively. Iter denotes the number of GMRES iterations and 
Time denotes the CPU time in seconds. GMRES, GSS, IGSS, F-GSS, F-IGSS, and mldivide 
denote GMRES with no preconditioning, GMRES preconditioned by the multistep GSS itera¬ 
tions, its inexact variant, FGMRES preconditioned by the multistep GSS iterations, its inexact 
variant, and the Matlab direct solver function mldivide, respectively. Hence, the CPU time 
for GMRES preconditioned by the multistep GSS iterations will improve with a sophisticated 
choice of the value of i. The least CPU time for each number of grids among the iterative 
methods is denoted by bold texts. 

The number of GSS and IGSS iterations was set to three throughout for simplicity, which 
is not necessarily optimal in terms of the CPU time. For example, GMRES preconditioned 
by six GSS iterations took 76.67 seconds to attain the stopping criterion for the problem with 
H = 10 -5 and q = 36. 

Table 4.1 shows that for well-conditioned problems /x = 1, IGSS (£ = 1) took the least CPU 
time to attain the stopping criterion among the iterative methods except for the small problem 
with grids 16 x 16. Table 4.2 shows that for ill-conditioned problems n = 10~ 5 , GSS (£ = 3) 
took the least CPU time among the iterative methods. For small problems with grids 16 x 16 
and 24 x 24, FGMRES took larger CPU time than other iterative methods. For ill-conditioned 
problems fx = 10~ 5 , although FGMRES required the fewest numbers of iterations, it did not 
outperform other methods in terms of the CPU time. The cost for solving the linear system 
Az = Vk with the stopping criterion ||u^ — Az^ +1 || < |cfc| (Theorem 3.1) was not marginal in 
FGMRES, since the value of |cfc| becomes small in the final FGMRES iterations. Note that the 
Matlab direct solver mldivide function gave more accurate solutions than the iterative methods 
within less CPU time. 


Table 4.1: Number of iterations and CPU time for (4.1) with fi = 1. 


Grids 

16 x 16 

24 x 24 

32 x 32 


Iter 

Time 

Iter 

Time 

Iter 

Time 

GMRES 

145 

0.279 

212 

0.321 

310 

2.505 

GSS (l = 1) 

19 

0.071 

20 

0.315 

23 

3.311 

GSS (£ = 3) 

13 

0.057 

15 

0.336 

17 

3.311 

IGSS (l = 1) 

18 

0.136 

19 

0.258 

21 

1.088 

IGGS {£ = 3) 

15 

0.170 

17 

0.699 

19 

2.966 

F-GSS 

29 

0.062 

37 

0.362 

40 

3.238 

F-IGSS 

29 

0.395 

38 

0.739 

39 

4.132 

mldivide 

0.011 

0.034 

0.066 


Table 4.2: Number of iterations and CPU time for (4.1) with /.i = 10 5 . 


Grids 

16 x 16 

24 x 24 

32 x 32 


Iter 

Time 

Iter 

Time 

Iter 

Time 

GMRES 

766 

2.915 

1,723 

17.51 

3,861 

391.3 

GSS {£ = 1) 

740 

2.693 

1,585 

16.32 

3,036 

275.2 

GSS {l = 3) 

561 

1.885 

1,130 

10.64 

1,549 

107.2 

IGSS (£ = 1) 

748 

5.668 

1,587 

42.62 

3,026 

406.6 

IGSS (£ = 3) 

594 

7.989 

1,155 

59.03 

1,586 

276.9 

F-GSS 

32 

17.17 

35 

60.29 

35 

451.0 

F-IGSS 

33 

295.5 

37 

500.8 

37 

4113. 

mldivide 

0.011 

0.015 

0.029 


4.2 Multistep Hermitian and skew-Hermitian splitting iteration precondi¬ 
tioning. 

We give numerical experiment results on singular and positive semidefinite linear system for 
( 1 . 1 ) with the generalized saddle point structure 


Ax 


' C 
- B J 



(4.3) 


where C € R pxp and G £ M. qxq are symmetric and positive semidefinite, and B £ MP xq . The 
results compare GMRES preconditioned by l iterations and FGMRES preconditioned by 1^ 
iterations of the Hermitian and skew-Hermitian splitting (HSS) [ 6 ] 


(cd + = (al — S)z W + v^, 

(ad + 5)z(* + h = (al — T^zh+V 2 ) + v^, 


(4.4) 


and its inexact variant (IHSS) with GMRES with no preconditioning and the standard HSS 
and IHSS preconditioning t = 1 for (4.3) and a sparse direct solver, where T~L = (A + AJ)/ 2 , 
S = (A — H t )/ 2, and a £ M. The former and latter systems of (4.4) were solved by using the 
Cholesky and LU factorizations, respectively, for HSS, and by using the conjugate gradient (CG) 
method [25] and the LSQR method [32], respectively, with the maximum number of iterations 
n, with the initial iterate equal to zero, and with the stopping criterion 10 -1 in terms of the 
relative residual 2-norm for the IHSS iterations [ ]. The maximum number of IHSS iterations 
for FGMRES was n. 

We generated test problems with the structure 


(U ffi V) J A{U 0 V) 


C 0 O 600 
-6 t 0O G0O 


b = A[l,2,...,n} J , 


(4.5) 
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where U € M pxp and V € R qxq are orthogonal matrices such that U T CU = C © O, V J GV = 
G © O, and U T BV = B © O. We set U T CU = diag(y?(l), y?(2),... ,<p(p — q — 1)) ® O € 
M pxp for <p(i) = k*/(p-9-!), j € N, U T GV = diag(^(l),^(2),... , ip(q - 2)) © O € R« x « for 
= K i ^ q ~ 1 '>, and V J B J U = [V T G T V,0] € R px? , where k = 10“U We set g = 16, 32, 
and 64, p = q 2 , nonzero density 0.1% of A, and j = 3, 6 , and 9 to show the effect of the 
condition number on the convergence. The value of j determines the condition number of A 
such as ||A||||At|| = \pl x 10 J , where A' is the pseudo inverse of A. The orthogonal matrices 
U 6 M pxp and V € R qxq were the products of random Givens rotations. Hence, the HSS iteration 
matrix (cel + 5) _1 (al — 7i)(al + — S) is semiconvergent [3, Theorem 3.6], and GMRES 

preconditioned by multistep HSS iterations determines a solution of (4.3) without breakdown 
(Theorem 2.4). On the other hand, the GMRES methods preconditioned by multistep IHSS 
and its iterations are not guaranteed to determine a solution without breakdown. 

Multistep HSS and IHSS iterations involve two parameters: the iteration parameter a and 
the number of HSS iterations l. Several techniques were proposed for estimating an optimal 
value of the HSS iteration parameter for the nonsingular case. As pointed out by a referee, 
parameter estimation techniques proposed in [5], [2], [26], [16] are developed for the present case. 
Huang’s technique need not modify for the singular case. In Chen’s technique, the minimum 
eigenvalue of % and the minimum singular value of S were replaced by the nonzero ones. After 
the value of the iteration parameter a was determined, the number of iterations i was determined 
by applying the HSS iterations alone to (4.3) and adopted the smallest between 10 and the 
smallest number of i which satisfies the relative difference norm ||z^ -1 ^ — || < 10 — 1 ||z®||. 

The CPU times required by Huang’s technique to determine the values of the HSS iteration 
parameter were 0.001 seconds for q = 16, 0.002 seconds for q = 32, and 0.019 seconds for q = 64. 
Bai et al.’s technique [5] and Bai’s technique [2] did not give more reasonable values of the HSS 
iteration parameter than Huang’s [26] and Chen’s [16] techniques for the test problems. Note 
that Bai’s technique [2] is for the saddle-point problem instead of the generalized saddle-point 
problem (4.3), and does not take into account the (2,2) block of (4.3) for the estimation. 

Tables 4.3-4.5 give the optimal and estimated values of the HSS iteration parameter and 
the value of the corresponding pseudo spectral radius of the HSS iteration matrix. The optimal 
value of the HSS iteration parameter a exp was experimentally determined to minimized the 
pseudo-spectral radius of the HSS iteration matrix. The values of the HSS iteration parameters 
which were estimated by using Huang’s and Chen’s techniques are denoted with subscript C and 
H, respectively. Chen’s technique estimated the values of the parameter close to the optimal 
one of the parameter which were experimentally determined. 

Tables 4.6-4 .8 give the number of the iterations and the CPU time in seconds for the test 
problems with different sizes and condition numbers. HSS, IHSS, HSS', IHSS', F-HSS', and F- 
IHSS' denote GMRES preconditioned by the HSS preconditioner, its inexact variant, GMRES 
preconditioned by the multistep HSS iterations, its inexact variants, FGMRES preconditioned 
by multistep HSS iterations, and its inexact variant, respectively, f means that CG or LSQR 
for the linear systems did not attain the stopping criterion within n iterations or the IHSS 
iterations did not satisfy ||u^ — A.z^|| < |cfc| within n iterations for the indicated number of 
iterations, f means that the Matlab direct solver mldivide function fails to give a solution, i.e., 
some of its entries are Not a Number (NaN). 

HSS' took the least CPU time to attain the stopping criterion among the iterative methods 
except for the case (j, q) = (9,16). Bai and Chen’s techniques tended to give reasonable values 
of the HSS iteration parameter for well-conditioned or small problems such as the cases (j, q) = 
(3,16), (3,64), (6,16), ( 6 , 32), whereas Huang’s technique tended to give reasonable values of the 
HSS iteration parameter for ill-conditioned or large problems such as the cases (j. q) = (6,64), 
(9,16), (9,32), (9,64). Although F-IHSS' took the fewest numbers of iterations, it did not 
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Table 4.3: HSS parameter a and pseudo spectral radius u(H(a)) of the iteration matrix of (4.4) 
for (4.3), j = 3. 


q 

16 

32 

64 

C^exp 

0.03162 

0.03162 

0.03162 

(t^exp)) 

0.93869 

0.93869 

0.93869 

oh 

0.15678 

0.08055 

0.03295 

u(H(a H )) 

0.98732 

0.97548 

0.94109 

ac 

0.03162 

0.03162 

0.03162 

v(H(a c )) 

0.93869 

0.93869 

0.93869 


Table 4.4: HSS parameter a and pseudo spectral radius v(H(a)) of the iteration matrix of (4.4) 
for (4.3), j = 6. 


q 

16 

32 

64 

C^exp 

0.00100 

0.00100 

0.00100 

(o; e xp)) 

0.99800 

0.99800 

0.99800 

a H 

0.14289 

0.08068 

0.03610 

v(H(a H )) 

0.99999 

0.99998 

0.99994 

ac 

0.00100 

0.00100 

0.00100 

v{H(a c )) 

0.99800 

0.99800 

0.99800 


Table 4.5: HSS parameter a and pseudo spectral radius v(H(a)) of (4.4) the iteration matrix 
for (4.3), j = 9. 


q 

16 

32 

64 

C^exp 

3.16e-5 

3.16e-5 

3.16e-5 

(oiexp)) 

0.99993 

0.99993 

0.99993 

a H 

0.14102 

0.13766 

0.03878 

v(H{an)) 

1.00000 

1.00000 

1.00000 

ac 

3.16e-5 

3.16e-5 

3.16e-5 

v(H(a c )) 

0.99993 

0.99993 

0.99993 


outperform other methods in terms of the CPU time. MSS' did not converge for all test 
problems. The Matlab direct solver mldivide function failed to give a solution for the large 
cases q = 64, although it outperformed the iterative methods for the other cases, except for the 
case j = 6. 

Comparing Tables 4.6-4.8 with Tables 4.3-4.5, we see that these estimated optimal values 
of the HSS iteration parameter in terms of the pseudo spectral radius did not give optimal CPU 
time for HSS'. This implies that a small pseudo-spectral radius does not necessarily gives a fast 
convergence of HSS, HSS', and IHSS (see also Theorem 2.5). 


5 Conclusions 

We considered applying several steps of matrix splitting iterations as a preconditioner to GM- 
RES and FGMRES for solving singular linear systems. We gave sufficient conditions such that 
GMRES and FGMRES preconditioned by multistep matrix splitting iterations determine a so¬ 
lution without breakdown, and a convergence bound of GMRES preconditioned by multistep 
matrix splitting iterations based on a spectral analysis. We presented a complexity issue of using 
multistep matrix splitting iteration preconditioning more than one step for GMRES. Numerical 
experiments showed that GMRES preconditioned by the multistep GSS and HSS iterations is 
efficient compared to previous methods including FGMRES for large and ill-conditioned prob- 
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Table 4.6: Number of iterations and CPU time for (4.3), (4.5) with j = 3 



<1 

16 

32 

64 



£ 

Iter 

Time 

£ 

Iter 

Time 

£ 

Iter 

Time 

GMRES 


112 

0.061 


159 

0.154 


167 

0.919 


HSS 


47 

0.018 


68 

0.049 


81 

0.493 


HSS' 

7 

19 

0.010 

10 

15 

0.039 

10 

15 

0.405 

C^exp 

IHSS 


52 

0.089 


56 

0.137 


70 

0.834 


IHSS' 

7 

29 

t 

10 

1 

t 

10 

1 

t 


F-HSS' 


19 

0.014 


1 

t 


1 

t 


F-IHSS' 


8 

0.267 


11 

0.349 


11 

1.442 


HSS 


63 

0.026 


62 

0.043 


79 

0.453 


HSS' 

5 

22 

0.011 

8 

15 

0.033 

10 

15 

0.405 

oh 

IHSS 


65 

0.083 


63 

0.129 


64 

0.808 

IHSS' 

5 

1 

t 

8 

1 

t 

10 

1 

t 


F-HSS' 


35 

0.017 


1 

t 


1 

t 


F-IHSS' 


14 

0.337 


15 

0.425 


13 

1.327 

mldivide 

0.000 

0.001 

fO.Oll 


Table 4.7: Number of iterations and CPU time for (4.3), (4.5) with j = 6 



<7 

16 

32 

64 



£ 

Iter 

Time 

£ 

Iter 

Time 

£ 

Iter 

Time 

GMRES 


163 

0.117 


464 

1.136 


1,014 

31.50 


HSS 


126 

0.073 


209 

0.262 


373 

4.634 


HSS' 

10 

49 

0.046 

10 

71 

0.164 

10 

109 

2.139 

C^exp 

IHSS 


130 

0.517 


262 

1.575 


738 

36.84 

etc 

IHSS' 

10 

1 

t 

10 

1 

t 

10 

1 

t 


F-HSS' 


1 

t 


1 

t 


1 

t 


F-IHSS' 


4 

t 


2 

t 


1 

t 


HSS 


162 

0.116 


436 

0.999 


326 

3.397 


HSS' 

5 

146 

0.136 

7 

150 

0.309 

10 

76 

1.520 

ttH 

IHSS 


163 

0.214 


444 

1.444 


362 

6.001 

IHSS' 

5 

1 

t 

7 

1 

t 

10 

1 

t 


F-HSS' 


32 

t 


1 

t 


1 

t 


F-IHSS' 


6 

t 


11 

t 


12 

t 

mldivide 

0.000 

$0,001 

$0,011 


Table 4.8: Number of iterations and CPU time for (4.3), (4.5) with j = 9 



Q 

16 

32 

64 



£ 

Iter 

Time 

£ 

Iter 

Time 

£ 

Iter 

Time 

GMRES 


117 

0.062 


353 

0.672 


852 

23.80 


HSS 


169 

0.125 


511 

1.333 


1,308 

49.83 


HSS' 

10 

162 

0.219 

10 

317 

1.063 

10 

487 

13.61 

<^exp 

IHSS 


107 

0.717 


297 

3.455 


938 

76.97 

etc 

IHSS' 

10 

1 

t 

10 

1 

t 

10 

1 

t 


F-HSS' 


1 

t 


1 

t 


1 

t 


F-IHSS' 


1 

t 


1 

t 


1 

t 


HSS 


117 

0.064 


348 

0.659 


346 

4.169 


HSS' 

5 

111 

0.089 

6 

190 

0.389 

10 

80 

1.503 

CtH 

IHSS 


117 

0.135 


350 

0.895 


381 

6.628 

IHSS' 

5 

1 

t 

6 

1 

t 

10 

1 

t 


F-HSS' 


30 

t 


1 

t 


1 

t 


F-IHSS' 


9 

t 


8 

t 


12 

t 

mldivide 

0.000 

0.001 

$0,011 
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Appendix 

If index(A) = minjd € No| rank A d = rank A d+l }, where A 0 = I and I is the identity matrix [12, 
Definition 7.2.1], then d > index(A) is equivalent to lZ(A d ) n Af(A d ) = {0} [12, Lemma 7.6.1]. 
The following theorem gives conditions such that GMRES determines a solution without break¬ 
down for 

index(A) > 1 . 

Theorem 5.1. GMRES determines a solution of Ax = b without breakdown for all b g lZ(A d ) 
and for all x$ € TZ(A d ~ 1 ) + AT (A) if and only if d> index(A). 

Proof. Assume d > index(A), or lZ(A d ) C\Af(A d ) = {0}. Let b € lZ(A d ) and xq € TZ(A d ~ 1 ) + 
Af{A). Then, r 0 € 7Z(A d ) and !C k C K{A d ). If k is the smallest positive integer such that 
dim/Cfc < k, then GMRES determines a solution of Ax = b at step k — 1 (see [11, Theorem 
2.2]). Now, assume dim/Q = i, 1 < i < k. Since /Q+i = r*oUA/Q, we have dim A 1C, = dim /C* = i 
for i = 1, 2 ,... , k — 1. Let the columns of V € W ixk form a basis of K, k - If dim AK, k < dim/Cfc, 
then there exists c / 0 such that AVc = 0. Since Vc / 0 for c / 0, we have KLkGAf(A) / {0}. 
From /Cfc C lZ(A d ) and Af(A) C Af(A d ), we have lZ(A d ) n A f{A d ) / {0}, which contradicts 
with d > index(A). Hence, dim AtC k = dim/C^ for all k g N. Since GMRES does not break 
down through rank deficiency of the least squares problem min zg ^; fe ||r*o — Az ||, the sufficiency 
is shown from [11, Theorem 2.2]. 

On the other hand, assume d < index(A). Then, A f{A d ) C A^A^" 1 " 1 ). There exits s / 0 such 
that s 0 Af(A d ) and s g A f{A d+1 ). Let t = A d s. Then, f ^ 0 and At = A d+l s = 0. Hence, 
there exists t / 0 such that t g lZ(A d ) 0 Af(A). Let b = t + Ax o for xq g 1Z(A d ~ l ) + Af(A). 
Then, b g 1Z(A d ) and = b — Ax o = t / 0. Since Aro = 0, we have rq = b — Ax\ = 
b — A(xq + cro) = tq — cAro = tq / 0 for c g R and dim A/Ci = 0 < dim/Ci = 1. Hence, 
GMRES breaks down at step 1 before determining a solution of Ax = b. Therefore, we complete 
the proof. □ □ 

This theorem agrees with [30, Theorem 2.2] for d = 1. Although similar results to Theorem 
5.1 were given in [12], no attention was paid there to the uniqueness of the GMRES iterate x k , 
i.e., the dimensions of AK. k and JC k . 
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